Große TODO Liste: gezeitenreibung_final.ipynb¶

gezeitenreibung_final.ipynb möglichst nicht ändern, sondern nur die Datei gezeitenreibung.ipynb bearbeiten. Nur fertige komplette Sachen in gezeitenreibung_final.ipynb kopieren.

TODOS:¶

--

Erst wenn oben alles fertig ist:¶

Zellen von gezeitenreibung.ipynb in gezeitenreibung_final.ipynb kopieren. $ \newline $ Baryzentrum in barycenter umbenenen, das wird aber im Augenblick zu verstreut verwendet.

$ \def\dv#1{ \frac{\mathrm{d}}{\mathrm{d} #1} } $ $ \def\dd#1{ \ \mathrm{d} #1 } $ $ \def\Vec#1{ \overrightarrow{#1} } $ $ \def\VecT#1#2{ \begin{pmatrix} #1 \\ #2 \end{pmatrix} } $ $ \def\VecD#1#2#3{ \begin{pmatrix} #1 \\ #2 \\ #3 \end{pmatrix} } $ $ \def\VecE#1{ \hat{e}_{#1} } $ $\providecommand{\e}[1]{\ensuremath{\cdot 10^{#1}}}$ $\providecommand{\fehlt}{\textcolor{red}{\textbf{Fehlt!\dots}}}$ $\providecommand{\todo}{\textcolor{red}{\textbf{\huge{ToDo}}}}$

In [ ]:
mErde, mOzean, RErde, mMond, rMondBahn, TMondBahn, G, TErdRotation, tau = 5.9721986*10**24, 0.0014*10**24, 6.3675*10**6, 7.3459*10**22, 3.836*10**8, 27.32166140*24*3600, 6.67430*10**(-11), 86164.100, 0.0021

import numpy as np
import matplotlib.pyplot as plt
from scipy.constants import pi
import scipy
from scipy.integrate import solve_ivp
from scipy.integrate._ivp.ivp import OdeResult # :) https://github.com/scipy/scipy/blob/main/scipy/integrate/_ivp/ivp.py
from matplotlib.animation import FuncAnimation
from IPython.display import Image #to display animations
from matplotlib import cm

# Beautiful plots
plt.rcParams['text.usetex'] = True
plt.rcParams['font.family'] = "serif"
plt.rcParams['font.serif'] = "Computer Modern Roman"

# Useful StackPosts
# https://stackoverflow.com/questions/59634279/solve-ivp-error-required-step-size-is-less-than-spacing-between-numbers

# Executing everything all the time takes too long
fastExecution = True
In [ ]:
# NICHT ANFASSEN! WENN, DANN UNBEDINGT IN FINAL AUCH ÄNDERN!
def baryzentrum(m1, m2, r): return r*m2/(m1+m2)

def iv_stable_orbit_2body():
    abstand_baryzentrum_erde = baryzentrum(mErde, mMond, rMondBahn)
    return [[-abstand_baryzentrum_erde, 0], [rMondBahn-abstand_baryzentrum_erde, 0], [0, -2*pi*abstand_baryzentrum_erde/TMondBahn], [0, 2*pi*(rMondBahn-abstand_baryzentrum_erde)/TMondBahn], [mErde, mMond]]

def eq_motion_2body(t, state, mass):
    x_1, y_1, x_2, y_2, vx_1, vy_1, vx_2, vy_2 = state
    dist_em = ((x_2 - x_1)**2 + (y_2 - y_1)**2 )**0.5
    ax_1 = G * mass[1] / (dist_em**3) * (x_2 - x_1)
    ay_1 = G * mass[1] / (dist_em**3) * (y_2 - y_1)
    ax_2 = G * mass[0] / (dist_em**3) * (x_1 - x_2)
    ay_2 = G * mass[0] / (dist_em**3) * (y_1 - y_2)
    return np.array([vx_1, vy_1, vx_2, vy_2, ax_1, ay_1, ax_2, ay_2])

def two_body_problem(pos_body_1: list, pos_body_2: list, vel_body_1: list, vel_body_2: list, mass: list, t_max: float, steps=1, method='RK45', t_start=0):
    solution = solve_ivp(eq_motion_2body, [t_start, t_max], [*pos_body_1, *pos_body_2, *vel_body_1, *vel_body_2], args=(mass,), atol=1e-4*rMondBahn, max_step=TMondBahn/1000, steps=steps, method=method)
    x1, y1, x2, y2, v_x1, v_y1, v_x2, v_y2 = solution.y
    return [solution.t, x1, y1, x2, y2]

def tide_acceleration(x_M, y_M, x_i, y_i, phi_i, ax_E, ay_E, omega_i, m_M):
    part_a, part_b = G*m_M* ((y_M - y_i)*np.cos(phi_i) - (x_M - x_i)*np.sin(phi_i)) / ((x_M - x_i)**2 + (y_M - y_i)**2)**(3/2), (RErde*omega_i**2 - (ax_E*np.cos(phi_i) + ay_E*np.sin(phi_i)))
    ax_i, ay_i = -part_a*np.sin(phi_i) - part_b*np.cos(phi_i), part_a*np.cos(phi_i) - part_b*np.sin(phi_i)
    return (-(ax_i - ax_E)*np.sin(phi_i) + (ay_i - ay_E)*np.cos(phi_i))/RErde

def eq_motion_4body(t, state, mass, calculate_tides=True):
    x_E, y_E, x_M, y_M, phi_1, phi_2, vx_E, vy_E, vx_M, vy_M, omega_1, omega_2 = state
    m_E, m_M, m_1, m_2 = mass
    x_1, y_1, x_2, y_2 = x_E + RErde*np.cos(phi_1), y_E + RErde*np.sin(phi_1), x_E + RErde*np.cos(phi_2), y_E + RErde*np.sin(phi_2)
    mu_x = m_E + m_1*np.cos(phi_1)**2 + m_2*np.cos(phi_2)**2
    mu_y = m_E + m_1*np.sin(phi_1)**2 + m_2*np.sin(phi_2)**2
    a = (G*m_M*m_E)/((x_M - x_E)**2 + (y_M - y_E)**2)**(3/2)
    b = m_1*np.sin(phi_1)*np.cos(phi_1) + m_2*np.sin(phi_2)*np.cos(phi_2)
    c1 = m_1*RErde*omega_1**2 + (G*m_M*m_1*((x_M - x_1)*np.cos(phi_1) + (y_M - y_1)*np.sin(phi_1)))/((x_M - x_1)**2 + (y_M - y_1)**2)**(3/2)
    c2 = m_2*RErde*omega_2**2 + (G*m_M*m_2*((x_M - x_2)*np.cos(phi_2) + (y_M - y_2)*np.sin(phi_2)))/((x_M - x_2)**2 + (y_M - y_2)**2)**(3/2)    
    ax_E = 1/(mu_x - b**2/mu_y) * (a*(x_M - x_E) -b/mu_y*(a*(y_M - y_E) + c1*np.sin(phi_1) + c2*np.sin(phi_2)) + c1*np.cos(phi_1) + c2*np.cos(phi_2))
    ay_E = 1/mu_y*(a*(y_M - y_E) - b*ax_E + c1*np.sin(phi_1) + c2*np.sin(phi_2))
    ax_M = G*((m_E*(x_E - x_M))/((x_E - x_M)**2 + (y_E - y_M)**2)**(3/2) + (m_1*(x_1 - x_M))/((x_1 - x_M)**2 + (y_1 - y_M)**2)**(3/2) + (m_2*(x_2 - x_M))/((x_2 - x_M)**2 + (y_2 - y_M)**2)**(3/2))
    ay_M = G*((m_E*(y_E - y_M))/((x_E - x_M)**2 + (y_E - y_M)**2)**(3/2) + (m_1*(y_1 - y_M))/((x_1 - x_M)**2 + (y_1 - y_M)**2)**(3/2) + (m_2*(y_2 - y_M))/((x_2 - x_M)**2 + (y_2 - y_M)**2)**(3/2))
    if not calculate_tides: return [vx_E, vy_E, vx_M, vy_M, omega_1, omega_2, ax_E, ay_E, ax_M, ay_M, 0, 0] # For code reuse with friction
    alpha_1, alpha_2 = tide_acceleration(x_M, y_M, x_1, y_1, phi_1, ax_E, ay_E, omega_1, m_M), tide_acceleration(x_M, y_M, x_2, y_2, phi_2, ax_E, ay_E, omega_2, m_M) 
    return[vx_E, vy_E, vx_M, vy_M, omega_1, omega_2, ax_E, ay_E, ax_M, ay_M, alpha_1, alpha_2]

def center_of_mass(mass, position_x , position_y):
    return [np.sum(mass*position_x)/np.sum(mass), np.sum(mass*position_y)/np.sum(mass)]

def iv_stable_orbit_4body():
    mass = [mErde, mMond, mOzean/2, mOzean/2] # Each of the tides gets half of the ocean's mass
    abstand_baryzentrum_erde = center_of_mass(mass, np.array([0, rMondBahn, -RErde, RErde]), np.array([0, 0, 0, 0]))[0]
    x0_Erde = [-abstand_baryzentrum_erde, 0]
    x0_Mond = [rMondBahn - abstand_baryzentrum_erde, 0]
    vy0_Erde = 2*pi*abstand_baryzentrum_erde/TMondBahn
    vy0_Mond = 2*pi*(rMondBahn-abstand_baryzentrum_erde)/TMondBahn
    v0_Erde = [0, -vy0_Erde]
    v0_Mond = [0, vy0_Mond]
    phi_1, phi_2 = 0, pi
    omega_1, omega_2 = 2*pi/TMondBahn, 2*pi/TMondBahn 
    return [x0_Erde, x0_Mond, [phi_1, phi_2], v0_Erde, v0_Mond, [omega_1, omega_2], mass]

def four_body_problem(pos_body_1: list, pos_body_2: list, phi_i: list, vel_body_1: list, vel_body_2: list, omega_i: list, mass: list, t_max: float, fun=eq_motion_4body, rtol=1e-12, atol=1e-12, t_start=0):
    solution = solve_ivp(fun=fun, t_span=[t_start, t_max], y0=[*pos_body_1, *pos_body_2, *phi_i, *vel_body_1, *vel_body_2, *omega_i], args=(mass,), rtol=rtol, atol=atol)
    x_E, y_E, x_M, y_M, phi_1, phi_2, vx_E, vy_E, vx_M, vy_M, omega_1, omega_2 = solution.y
    return [solution.t, x_E, y_E, x_M, y_M, phi_1, phi_2, omega_1, omega_2]

Zweikörperproblem¶

In [ ]:
# Calculate the solution of the two body problem
t, x_1, y_1, x_2, y_2 = two_body_problem(*iv_stable_orbit_2body(), 3*TMondBahn)
/home/finn/.local/lib/python3.10/site-packages/scipy/integrate/_ivp/common.py:39: UserWarning: The following arguments have no effect for a chosen solver: `steps`.
  warn("The following arguments have no effect for a chosen solver: {}."
  • Vergleichen Sie alle Ergebnisse auch mit der exakten Lösung des Zweikörperproblems.
Vergleich mit der exakten Lösung des Zweikörperproblems:¶

Eine exakte analytische Lösung des Zweikörperproblems ist mit genügend Annahmen möglich. Dazu wird im Schwerpunktssytem gearbeitet. Aus den Bewegungsgliechungen der Körper die nur mit der Schwerkraft wechselwirken, kann man folgern, dass der Schwerpunkt des Systems eine gleichförmige Bewegung ausführt. Mit unseren Anfangsbedingungen bewegt dieser sich aber nicht.

Die beiden Körper können sich um ihren Schwerpunkt je nach Anfangsbedingungen auf verschiedenen Bahnen bewegen. Dies lässt sich durch die Exzentrizität beschreiben. Mögliche Bahnen sind Kreise, Ellipsen, Parabeln und Hyperbeln. Wir haben die Anfangsbedingungen so gewählt, dass die Bahn eine Kreisbahn ist. $\newline$ Unsere Ergebnise stimmen also mit der exakten Lösung des Zweikörperproblems überein, da wir zwei geschlossene Kreisbahnen um den gemeinsamen Schwerpunkt (mit unseren Anfangsbedingungen der Ursprung) erhalten. Diese Kreisbahnen lassen sich mit einfachen trigonometrischen Funktionen beschreiben. \begin{align} r_E(t) = \VecT{x_E(t)}{y_E(t)} = \VecT{r_{E_0}\cos(\omega_E t)}{r_{E_0}\sin(\omega_E t)} \\ r_M(t) = \VecT{x_M(t)}{y_M(t)} = \VecT{r_{M_0}\cos(\omega_M t)}{r_{M_0}\sin(\omega_M t)} \end{align} Diese lassen sich ebenfalls leicht plotten.

Unser numerischer Ansatz ist natürlich in der Lage alle möglichen Bahnen zu berechnen, also auch Ellipsen und Hyperbeln, es ließe sich also auch ein Komet oder ein Sling-Shot-Manöver berechnen. Mit den geeigneten Anfangsbedingungen ließen sich diese auch analytisch berechnen.

Eine weitere Möglichkeit unsere Ergebnisse zu überprüfen ist das dritte Kepplersche Gesetz, es besagt, dass die Umlaufzeit $T$ des Zweikörpersystem beträgt: $$T^2 = \frac{4\pi^2 \cdot a^3}{G\cdot M} $$ Mit $G$ als Gravitationskonstante, $M = m_{Erde}+m_{Mond}$ und $a$ als großer Halbachse welche in unserem genäherten Fall der Kreisbahn $a=r_{MondBahn}$ entspricht.

In [ ]:
r_E0, r_M0 = iv_stable_orbit_2body()[0][0], iv_stable_orbit_2body()[1][0]
time_points = np.linspace(0, 3*TMondBahn, 1000)
x_E, y_E = np.vectorize(lambda t: r_E0 * np.cos(2*pi/TMondBahn * t))(time_points), np.vectorize(lambda t: r_E0 * np.sin(2*pi/TMondBahn * t))(time_points)
x_M, y_M = np.vectorize(lambda t: r_M0 * np.cos(2*pi/TMondBahn * t))(time_points), np.vectorize(lambda t: r_M0 * np.sin(2*pi/TMondBahn * t))(time_points)

# Plot the analytical solution with our boundary conditions
plt.rc('axes', axisbelow=True)
plt.figure(figsize=(5, 5))
plt.title("Analytische Bahnen des 2-Körper-Systems")
plt.grid('minor', 'minor', linestyle='-', linewidth=0.2)
plt.grid('major', 'major', linestyle='-', linewidth=0.8)
plt.minorticks_on()
plt.xlabel("x")
plt.ylabel("y")
plt.plot(x_E, y_E, 'b', linewidth=1, label='Erde')
plt.plot(x_M, y_M, 'r', linewidth=1, label='Mond')
plt.legend()
plt.show()

print(f'3. Keplersches Gesetze T: {TMondBahn} = {(4*pi**2*rMondBahn**3/(G*(mErde+mMond)))**(1/2)}')
3. Keplersches Gesetze T: 2360591.54496 = 2350028.5085322363

Vierkörperproblem mit Randbedingungen, intrinsischer Rotation und Reibung¶

Im Folgenden wollen wir das 4-Körper-Problem um die Reibung des Ozeans (der Flutberge) mit der Erde erweitern. Dies verlangsamt die intrinsische Erdrotation und sorgt für länger werdende Tage und eine größere Umlaufbahn des Mondes. Um dies zu beschreiben müssen wir mehrere zusätzliche Komponenten herleiten und mit ihnen unser Modell erweitern.

Die Reibung beschreiben wir zunächst mit der Reibungskonstante $ k = 2 \cdot 10^{-12}\,\frac{1}{\mathrm m}$.

Herleitung der zusätzlichen Bewegungsgleichung und der Reibungskraft¶

Wir wollen als erstes die Reibungskraft zwischen der Erde und den Flutbergen bestimmen. Wir wählen dafür folgenden Ansatz: $ \vec F_{\mathrm R,i}(\vec v_i)=-km_i|\vec v_i| \vec v_i $. Hierbei ist $\vec v_i$ die Geschwindigkeit des $i$-ten Flutbergs relativ zur rotierenden Erdoberfläche und $k$ eine effektive Reibungskonstante.

Die Geschwindigkeit des Flutbergs bestimmen wir über seine Winkelgeschwindikeit relativ zur Erdoberfläche $(\omega_E - \omega_i)$, mit der er sich um die z-Achse dreht. Hierbei ist $\omega_E $ die Winkelgeschwindigkeit der Eigenrotation der Erde und $\omega_i$ die Winkelgeschwindigkeit des Flutbergs.

Die Geschwindigkeit ist allgemein definiert als $\vec v_i = \vec \omega_i \times \vec r_i$: \begin{equation*} \vec v_i = (\omega_E - \omega_i) \cdot \VecD{0}{0}{1} \times r_{E} \cdot \vec e_{r, i} = \VecD{0}{0}{\omega_E - \omega_i} \times \VecD{r_{E} \cdot \cos(\phi_i)}{r_{E} \cdot \sin(\phi_i)}{0} = \VecD{-r_{E} \cdot \sin(\phi_i) \cdot (\omega_E - \omega_i)}{r_{E} \cdot \cos(\phi_i) \cdot (\omega_E - \omega_i)}{0} \end{equation*} Damit berechnen wir den Betrag: \begin{equation*} |\vec v_i| = \sqrt{r_{E}^2 \cdot (\omega_E - \omega_i)^2 \cdot \sin^2(\phi_i) + r_{E}^2 \cdot (\omega_E - \omega_i)^2 \cdot \cos^2(\phi_i)} = r_{E} \cdot |\omega_E - \omega_i| \end{equation*} Also ist die Reibungskraft: \begin{equation*} \vec F_{\mathrm R,i}(\vec v_i) = -km_i |\vec v_i| \vec v_i = -km_i r_{E} \cdot |\omega_E - \omega_i| \cdot \VecD{-r_{E} \cdot \sin(\phi_i) \cdot (\omega_E - \omega_i)}{r_{E} \cdot \cos(\phi_i) \cdot (\omega_E - \omega_i)}{0} = -km_i r_{E}^2 \cdot |\omega_E - \omega_i|(\omega_E - \omega_i) \cdot \VecD{- \sin(\phi_i)}{\cos(\phi_i)}{0} \end{equation*} Durch teilen durch die Masse des i.ten Flutbergs erhalten wir die Beschleuningung, die wir auf die Differentialgleichung des i.ten Flutbergs aus dem 4-Körper-Problem addieren.

Als nächstes wollen wir die Differentialgleichung für die intrinsische Rotation der Erde herleiten. Es gilt die Drehimpulsrelation $L = I \cdot \omega$. Wir nehmen die Erde als Kugel mit homogener Massenverteilung an, damit hat sie ein Trägheitsmoment von $I = \frac{2}{5} m_E r_E^2$. Für Drehimpuls und Drehmoment folgt: \begin{gather*} \vec L = \frac{2}{5} m_E r_E^2 \vec \omega_E \\ \vec M = \frac{\mathrm d \vec L}{\mathrm d t} = \frac{2}{5} m_E r_E^2 \vec{\dot \omega_E} \end{gather*} Das Drehmoment welches durch die Reibungskräfte entsteht ist: \begin{equation*} \vec M_{\mathrm R} = \sum_{i=1}^2 \vec r_i \times \vec F_{\mathrm R,i} \end{equation*} Dies berechnet sich für zwei Flutberge zu: \begin{equation*} \vec M_{\mathrm R} = \VecD{0}{0}{-k r_E^3 \left(m_1 \cdot |\omega_E - \omega_1|(\omega_E - \omega_1) + m_2 \cdot |\omega_E - \omega_2|(\omega_E - \omega_2) \right)} \end{equation*} Setzen wir die z-Komponenten der Drehmomente gleich (nur diese sind $\neq 0$), erhalten wir folgende Bewegungsgleichung für die intrinsische Rotation der Erde: \begin{gather*} \frac{2}{5} m_E r_E^2 \dot \omega_E = -k r_E^3 \left(m_1 \cdot |\omega_E - \omega_1|(\omega_E - \omega_1) + m_2 \cdot |\omega_E - \omega_2|(\omega_E - \omega_2) \right) \\ \dot \omega_E = -\frac{5 k r_E}{2 m_E} \left(m_1 \cdot |\omega_E - \omega_1|(\omega_E - \omega_1) + m_2 \cdot |\omega_E - \omega_2|(\omega_E - \omega_2) \right) \end{gather*}

Um dies im Code umzusetzen verwenden wir, wenn es geht, die schon implementierten Funktionen. Daher wird hier der selbe Ansatz wie zuvor gewählt, nur die Beschleunigung auf die Flutberge wird wie oben erklärt erweitert, und die zusätzliche Differentialgleichung wird hinzugefügt.

In [ ]:
def iv_stable_orbit_4body_friction():
    '''Initial conditions for the Earth-Moon-Tides system in a stable orbit'''
    x0_Erde, x0_Mond, phi_i, v0_Erde, v0_Mond, omega_i, mass = iv_stable_orbit_4body()
    phi0_E = 0 # The starting point of the earth's rotation is unconsequentially defined as 0
    omega0_E = 2*pi / TErdRotation # Angular velocity of the earth
    return [x0_Erde, x0_Mond, phi_i, phi0_E, v0_Erde, v0_Mond, omega_i, omega0_E, mass]

def tides_acceleration_friction(x_M, y_M, x_i, y_i, phi_i, ax_E, ay_E, omega_i, omega_E, m_M, k):
    '''Calculates the angular acceleration of a tide, adding the frictional force caused by the earth's rotation'''
    ax_friction = k * RErde**2 * np.abs(omega_i - omega_E)* (omega_i - omega_E) * np.sin(phi_i)
    ay_friction = -k * RErde**2 * np.abs(omega_i - omega_E)* (omega_i - omega_E) * np.cos(phi_i)
    part_a = G*m_M* ((y_M - y_i)*np.cos(phi_i) - (x_M - x_i)*np.sin(phi_i)) / ((x_M - x_i)**2 + (y_M - y_i)**2)**(3/2)
    part_b = (RErde*omega_i**2 - (ax_E*np.cos(phi_i) + ay_E*np.sin(phi_i)))
    ax_i = -part_a*np.sin(phi_i) - part_b*np.cos(phi_i) + ax_friction
    ay_i =  part_a*np.cos(phi_i) - part_b*np.sin(phi_i) + ay_friction

    alpha_i = (-(ax_i - ax_E)*np.sin(phi_i) + (ay_i - ay_E)*np.cos(phi_i))/RErde
    return alpha_i

def eq_motion_4body_friction(t, state, mass, k):
    x_E, y_E, x_M, y_M, phi_1, phi_2, phi_E, vx_E, vy_E, vx_M, vy_M, omega_1, omega_2, omega_E = state
    x_1, y_1, x_2, y_2 = x_E + RErde*np.cos(phi_1), y_E + RErde*np.sin(phi_1), x_E + RErde*np.cos(phi_2), y_E + RErde*np.sin(phi_2)
    m_E, m_M, m_1, m_2 = mass
    
    # Reuse the equations of motion from the 4-body problem for the earth and the moon:
    _, _, _, _, _, _, ax_E, ay_E, ax_M, ay_M, _, _ = eq_motion_4body(t, np.array([x_E, y_E, x_M, y_M, phi_1, phi_2, vx_E, vy_E, vx_M, vy_M, omega_1, omega_2]), mass, calculate_tides=False)

    # Tides angular acceleration with friction:
    alpha_1 = tides_acceleration_friction(x_M, y_M, x_1, y_1, phi_1, ax_E, ay_E, omega_1, omega_E, m_M, k)
    alpha_2 = tides_acceleration_friction(x_M, y_M, x_2, y_2, phi_2, ax_E, ay_E, omega_2, omega_E, m_M, k)

    # The earth's angular acceleration:
    alpha_E = (5*k*RErde)/(2*m_E) * (m_1*np.abs(omega_1 - omega_E)*(omega_1 - omega_E) + m_2*np.abs(omega_2 - omega_E)*(omega_2 - omega_E))

    return[vx_E, vy_E, vx_M, vy_M, omega_1, omega_2, omega_E, ax_E, ay_E, ax_M, ay_M, alpha_1, alpha_2, alpha_E]

def four_body_problem_friction(pos_body_1: list, pos_body_2: list, phi_i: list, phi_E: float, vel_body_1: list, vel_body_2: list, omega_i: list, omega_E: float, mass: list, k: float, t_max: float, rtol=1e-6, atol=1e-6, t_start=0):
    solution = solve_ivp(fun=eq_motion_4body_friction, t_span=[t_start, t_max], y0=[*pos_body_1, *pos_body_2, *phi_i, phi_E, *vel_body_1, *vel_body_2, *omega_i, omega_E], args=(mass, k), rtol=rtol, atol=atol)
    x_E, y_E, x_M, y_M, phi_1, phi_2, phi_E, vx_E, vy_E, vx_M, vy_M, omega_1, omega_2, omega_E = solution.y
    return [solution.t, x_E, y_E, x_M, y_M, phi_1, phi_2, phi_E, omega_1, omega_2, omega_E]

Graphische Darstellung und Vergleich zu ohne Reibung¶

In [ ]:
if not fastExecution:
    # 4 body without friction: 
    t, x_E, y_E, x_M, y_M, phi_1, phi_2, _, _ = four_body_problem(*iv_stable_orbit_4body(), 0.75*TMondBahn, rtol=1e-12, atol=1e-12) # 0.3s
    x_1, y_1, x_2, y_2 = x_E + RErde*np.cos(phi_1), y_E + RErde*np.sin(phi_1), x_E + RErde*np.cos(phi_2), y_E + RErde*np.sin(phi_2)
    t_l, x_E_l, y_E_l, x_M_l, y_M_l, phi_1_l, phi_2_l, omega_1_l, omega_2_l = four_body_problem(*iv_stable_orbit_4body(), 100*365*24*3600, rtol=1e-9, atol=1e-9)
    x_1_l, y_1_l, x_2_l, y_2_l = x_E_l + RErde*np.cos(phi_1_l), y_E_l + RErde*np.sin(phi_1_l), x_E_l + RErde*np.cos(phi_2_l), y_E_l + RErde*np.sin(phi_2_l)
    # 4 body with friction short time: 0.3s
    t_f, x_E_f, y_E_f, x_M_f, y_M_f, phi_1_f, phi_2_f, phi_E_f, omega_1_f, omega_2_f, omega_E_f = four_body_problem_friction(*iv_stable_orbit_4body_friction(), k=2e-12, t_max=0.75*TMondBahn, rtol=1e-12, atol=1e-12)
    x_1_f, y_1_f, x_2_f, y_2_f = x_E_f + RErde*np.cos(phi_1_f), y_E_f + RErde*np.sin(phi_1_f), x_E_f + RErde*np.cos(phi_2_f), y_E_f + RErde*np.sin(phi_2_f)
    # 4 body with friction  for 100 years: 11min
    t_fl, x_E_fl, y_E_fl, x_M_fl, y_M_fl, phi_1_fl, phi_2_fl, phi_E_fl, omega_1_fl, omega_2_fl, omega_E_fl = four_body_problem_friction(*iv_stable_orbit_4body_friction(), k=2e-12, t_max=100*365*24*3600, rtol=1e-12, atol=1e-12)
    x_1_fl, y_1_fl, x_2_fl, y_2_fl = x_E_fl + RErde*np.cos(phi_1_fl), y_E_fl + RErde*np.sin(phi_1_fl), x_E_fl + RErde*np.cos(phi_2_fl), y_E_fl + RErde*np.sin(phi_2_fl)
In [ ]:
if not fastExecution:
    # Plot the solution of the four body problem as a phase space diagram
    fig, axes = plt.subplots(nrows=5, ncols=4, figsize=(12, 12)) # Create a figure and a set of subplots
    fig.suptitle("Phasenraumdiagramme des Erde-Mond-Flutberg-Systems: ohne Reibung und mit Reibung")

    def init_phase_space(ax, a, b, title, xlabel, ylabel, color):
        ax.plot(a, b, color=color)
        ax.set_xlabel(xlabel)
        ax.set_ylabel(ylabel)
        ax.set_title(title)

    init_phase_space(axes[0,0], t, x_E, "x-Position Erde", "$t$ / s", "$x_1$ / m", "blue")
    init_phase_space(axes[0,1], t, x_M, "x-Position Mond", "$t$ / s", "$x_2$ / m", "red")
    init_phase_space(axes[1,0], t, y_E, "y-Position Erde", "$t$ / s", "$y_1$ / m", "blue")
    init_phase_space(axes[1,1], t, y_M, "y-Position Mond", "$t$ / s", "$y_2$ / m", "red")
    init_phase_space(axes[2,0], x_E, y_E, "Rotation Erde", "$x_E$ / m", "$y_E$ / m", "blue")
    init_phase_space(axes[2,1], x_M, y_M, "Rotation Mond", "$x_M$ / m", "$y_M$ / m", "red")
    init_phase_space(axes[3,0], t, phi_1, "Winkel Flutberg 1", "$t$ / s", "$\phi_1$ / rad", "blue")
    init_phase_space(axes[3,1], t, phi_2, "Winkel Flutberg 2", "$t$ / s", "$\phi_2$ / rad", "red")
    init_phase_space(axes[4,0], x_1, y_1, "Rotation Flutberg 1", "$x_1$ / m", "$x_2$ / m", "blue")
    init_phase_space(axes[4,1], x_2, y_2, "Rotation Flutberg 2", "$y_1$ / m", "$y_2$ / m", "red")

    init_phase_space(axes[0,2], t_f, x_E_f, "x-Position Erde mit Reibung", "$t$ / s", "$x_1$ / m", "blue")
    init_phase_space(axes[0,3], t_f, x_M_f, "x-Position Mond mit Reibung", "$t$ / s", "$x_2$ / m", "red")
    init_phase_space(axes[1,2], t_f, y_E_f, "y-Position Erde mit Reibung", "$t$ / s", "$y_1$ / m", "blue")
    init_phase_space(axes[1,3], t_f, y_M_f, "y-Position Mond mit Reibung", "$t$ / s", "$y_2$ / m", "red")
    init_phase_space(axes[2,2], x_E_f, y_E_f, "Rotation Erde mit Reibung", "$x_E$ / m", "$y_E$ / m", "blue")
    init_phase_space(axes[2,3], x_M_f, y_M_f, "Rotation Mond mit Reibung", "$x_M$ / m", "$y_M$ / m", "red")
    init_phase_space(axes[3,2], t_f, phi_1_f, "Winkel Flutberg 1 mit Reibung", "$t$ / s", "$\phi_1$ / rad", "blue")
    init_phase_space(axes[3,3], t_f, phi_2_f, "Winkel Flutberg 2 mit Reibung", "$t$ / s", "$\phi_2$ / rad", "red")
    init_phase_space(axes[4,2], x_1_f, y_1_f, "Rotation Flutberg 1 mit Reibung", "$x_1$ / m", "$x_2$ / m", "blue")
    init_phase_space(axes[4,3], x_2_f, y_2_f, "Rotation Flutberg 2 mit Reibung", "$y_1$ / m", "$y_2$ / m", "red")

    plt.tight_layout()
    plt.savefig("4body_friction_comp0.png", dpi=300)
    plt.close()

    fig, axes = plt.subplots(nrows=2, ncols=3, figsize=(15, 10))
    plt.rc('axes', axisbelow=True)
    plt.suptitle("Verlgeich 4 Körper Problem mit Reibung / ohne Reibung")

    def plot_subplot(ax, lines, title, xlabel, ylabel):
        for x, y, label, color in lines:
            ax.plot(x, y, color, linewidth=1, label=label)
        ax.set_xlabel(xlabel)
        ax.set_ylabel(ylabel)
        ax.set_title(title)

        ax.grid('minor', 'minor', linestyle='-', linewidth=0.2)
        ax.grid('major', 'major', linestyle='-', linewidth=0.8)
        ax.minorticks_on()
        ax.legend()

    # Distance from the origin for earth and moon
    r_Earth = np.sqrt(x_E_l**2 + y_E_l**2)
    r_Moon = np.sqrt(x_M_l**2 + y_M_l**2)
    r_Earth_f = np.sqrt(x_E_fl**2 + y_E_fl**2)
    r_Moon_f = np.sqrt(x_M_fl**2 + y_M_fl**2)

    plot_subplot(axes[0][0], [(x_E_l*10, y_E_l*10, 'Erde (*10)', 'k'), (x_M_l, y_M_l, 'Mond', 'r'), (x_1_l*10, y_1_l*10, 'Flutberg 1 (*10)', 'b'), (x_2_l*10, y_2_l*10, 'Flutberg 2 (*10)', 'b')], "Bahnen des 4-Körper Problems ohne Reibung", "x / m", "y / m")
    plot_subplot(axes[0][1], [(t_l, r_Moon, 'Mond', 'r')], "Abstand des Mondes zum Ursprung ohne Reibung", "t / s", "r / m")
    plot_subplot(axes[0][2], [(t_l, r_Earth, 'Erde', 'k')], "Abstand der Erde zum Ursprung ohne Reibung", "t / s", "r / m")
    plot_subplot(axes[1][0], [(x_E_fl*10, y_E_fl*10, 'Erde (*10)', 'k'), (x_M_fl, y_M_fl, 'Mond', 'r'), (x_1_fl*10, y_1_fl*10, 'Flutberg 1 (*10)', 'b'), (x_2_fl*10, y_2_fl*10, 'Flutberg 2 (*10)', 'b')], "Bahnen des 4-Körper Problems mit Reibung", "x / m", "y / m")
    plot_subplot(axes[1][1], [(t_fl, r_Moon_f, 'Mond', 'r')], "Abstand des Mondes zum Ursprung mit Reibung", "t / s", "r / m")
    plot_subplot(axes[1][2], [(t_fl, r_Earth_f, 'Erde', 'k')], "Abstand der Erde zum Ursprung mit Reibung", "t / s", "r / m")
    plt.tight_layout()
    plt.savefig("4body_friction_comp1.png", dpi=300)
    plt.close()

    fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(10, 5))
    plt.rc('axes', axisbelow=True)
    plt.suptitle("Verlgeich 4 Körper Problem mit Reibung / ohne Reibung")

    plot_subplot(axes[0][0], [(t_l, omega_1_l, '$\omega_1$ Flutberg 1', 'k'), (t_l, omega_2_l, '$\omega_2$ Flutberg 2', 'r'), (t_l, [2*pi/TMondBahn]*np.shape(t_l)[0], 'Startwinkelgeschwindigkeit Flutberg', 'g')], "Winkelgeschwindigkeiten der Flutberge ohne Reibung", "t / s", "$\omega$ / rad/s")
    plot_subplot(axes[0][1], [(t_l, [2*pi/TErdRotation]*np.shape(t_l)[0], '$\omega_E$ Winkelgeschwindigkeit Erde', 'k')], "Winkelgeschwindigkeit der Erde ohne Reibung", "t / s", "$\omega$ / rad/s")
    plot_subplot(axes[1][0], [(t_fl, omega_1_fl, '$\omega_1$ Flutberg 1', 'k'), (t_fl, omega_2_fl, '$\omega_2$ Flutberg 2', 'r'), (t_fl, [2*pi/TMondBahn]*np.shape(t_fl)[0], 'Startwinkelgeschwindigkeit Flutberg', 'g')], "Winkelgeschwindigkeiten der Flutberge mit Reibung", "t / s", "$\omega$ / rad/s")
    plot_subplot(axes[1][1], [(t_fl, omega_E_fl, '$\omega_E$ Winkelgeschwindigkeit Erde', 'k'), (t_fl, [2*pi/TErdRotation]*np.shape(t_fl)[0], 'Startwinkelgeschwindigkeit Erde', 'g')], "Winkelgeschwindigkeit der Erde mit Reibung", "t / s", "$\omega$ / rad/s")

    plt.tight_layout()
    plt.savefig("4body_friction_comp2.png", dpi=300)
    plt.close()

# Display image set the size to be 500x500 pixels
display(Image(filename='4body_friction_comp0.png', width=1000, height=2000))
display(Image(filename='4body_friction_comp1.png', width=1000, height=500))
display(Image(filename='4body_friction_comp2.png', width=1000, height=500))

Im Fall ohne Reibung ist $\omega_\mathrm{E}$ natürlich konstant, und im Fall mit Reibung nimmt $\omega_\mathrm{E}$ mit der Zeit ab, das heißt, dass die Tage auf der Erde länger werden. Hier sieht die Abnahme linear aus, es handelt sich aber um einen exponentiellen Abfall, welcher aufgrund der betrachteten Zeitskala aber linear erscheint.
Zudem ist zu sehen, dass Erde und Mond sich mit Reibung von ihrem gemeinsamen Schwerpunkt wegbewegen, was ohne Reibung nicht passiert. Das ist eine Folge davon, dass die Erdrotation schneller als die Umlaufbahn des Mondes ist. Die Flutberge werden aufgrund der Reibung vor die Verbindungslinie zwischen Erde und Mond gedrückt, was die Wirkung der Flutberge auf den Mond ändert. Der Mond erfährt eine Radialbeschleunigung, wodurch er auf eine weiter entfernte Bahn gelangt.
Vergleicht man die Winkelgeschwindigkeiten der Flutberge, so erkennt man, dass diese in beiden Fällen mit einer schwach modulierten Oszillation starten. Mit der Zeit schwächt sich dies mit Reibung ab, ohne Reibung bleibt der Effekt zeitlich gleich. Die Bewegung ist Folge davon, dass die Umlaufbahn nicht perfekt ein Kreis ist (Wie man auch in den Plots vom Abstand zum Uspung erkennen kann), weshalb der Mond die Flutberge zu unterschiedlichen Zeitpunkten unterschiedlich stark beschleunigt. Mit Reibung wird diese Bewegung abgebremst (Energiedissipation), weshalb die Amplitude der Oszillation abnimmt. Die Amplitude ist zu Beginn jedoch größer, weil die Flutberge durch die Reibung in Richtung der Erdrotation beschleunigt werden.

Wir stellen noch dar, wie sich der Mond von der Erde entfernt:

In [ ]:
#Show the incerase of the moons orbit
if not fastExecution:
    #Moon as seen from Earth:
    x = x_M_fl - x_E_fl
    y = y_M_fl - y_E_fl

    fig = plt.figure(figsize=(12,6))
    fig.suptitle('Vergrößerung der Mondbahn über 100 Jahre')
    ax = fig.add_subplot(121)
    ax2 = fig.add_subplot(122)

    ax.set_title('Mond im Erdsystem')
    ax.set_xlabel("$x$ / m")
    ax.set_ylabel("$y$ / m")
    ax.set_aspect('equal')
    ax.plot(x, y, color='r', linewidth=0.5, label='Mond')

    rect = plt.Rectangle((3.75e8, -1e5), 1.5e7, 1e7, facecolor="black", alpha=0.7, label='Ausschnitt')
    ax.add_patch(rect)
    ax.legend()

    ax2.set_title('Ausschnitt')
    ax2.set_xlabel("$x$ / m")
    ax2.set_ylabel("$y$ / m")
    ax2.set_aspect('equal')
    ax2.plot(x, y, 'r', linewidth=0.1, label='Mond')
    ax2.scatter(x[0], y[0], c='k', s=10, label='Startpunkt')
    ax2.set_xlim([3.8358e8, 3.8391e8])
    ax2.set_ylim([-1e5, 1e5])
    ax2.legend()

    plt.tight_layout()
    plt.savefig("4body_friction_Moon.png", dpi=300)
    plt.close()

display(Image(filename='4body_friction_Moon.png', width=1000, height=2000))

Berechnung der Tageslängenänderung¶

Für die erstmalige Berechnung der Tageslängenänderung verwenden wir die Reibungskonstante $k = 2 \cdot 10^{-12}\,\frac{1}{\mathrm m}$ und Massen für Flutberge von $ m_1 = m_2 = \frac{1}{2} m_{\mathrm{Ozean}} = \frac{1}{2} \cdot 0.0014 \cdot 10^{24} $

Die Verlangsamung der Erdrotation kann man auch in Realität beobachten. Deshalb vergleichen wir unsere Ergebnisse mit dem Literaturwert für die Tageslängenänderung über 100 Jahre $ \tau = 0.0021 \frac{\mathrm s}{\mathrm{100 a}}$.

Wie im Graphen für die Winkelgeschwindigkeiten der Flutberge zu erkennen ist braucht das System eine gewisse Einschwingzeit. Wir wählen konservativ 200 Jahre als Einschwingzeit (mit Hilfe des Graphen bestimmt). Im folgenden berechnen wir die Tageslängenänderung für 500 Jahre mit 100 Jahr Zeitschritten.

Die Simulation liefert uns die Winkelgeschwindigkeit der Erde $\omega_E$ zu gewissen Zeitpunkten. $ \newline $ Diese kann mit der Formel $\omega_E = \frac{2 \pi}{T} \Rightarrow T = 2 \pi \omega$ in die Tageslänge $T$ umgerechnet werden. Die Differenz der Tageslängen ergibt die Tageslängenänderung $ \tau $.

In [ ]:
if not fastExecution:
    year_in_seconds = 365*24*3600
    t_f5ha, _, _, _, _, _, _, _, _, _, omega_E_f5ha = four_body_problem_friction(*iv_stable_orbit_4body_friction(), k=2e-12, t_max=500*year_in_seconds, rtol=1e-12, atol=1e-12)
    angular_velocities = [omega_E_f5ha[i] for i in [np.argmin(np.abs(t_f5ha - j*year_in_seconds)) for j in [0, 100, 200, 300, 400, 500]]]
    tau_calc = [2*pi/angular_velocities[i+1] - 2*pi/angular_velocities[i] for i in range(5)]
    for i in range (5):
        print(f'τ für ({i}ha-{i+1}ha): {tau_calc[i]:.3f} s/ha')
else:
    print("τ für (0ha-1ha): 137.090 s/ha\nτ für (1ha-2ha): 137.086 s/ha\nτ für (2ha-3ha): 137.082 s/ha\nτ für (3ha-4ha): 137.078 s/ha\nτ für (4ha-5ha): 137.073 s/ha")
τ für (0ha-1ha): 137.090 s/ha
τ für (1ha-2ha): 137.086 s/ha
τ für (2ha-3ha): 137.082 s/ha
τ für (3ha-4ha): 137.078 s/ha
τ für (4ha-5ha): 137.073 s/ha

Das berechnete Ergebnis von ca. $ 137 \frac{\mathrm s}{\mathrm{100 a}}$ weicht stark vom Literaturwert $\tau$ ab.

Es ist deutlich zu erkennen, dass die Tageslängenänderung mit der Zeit abnimmt, da die Reibungskraft mit der Zeit abnimmt. Auf unseren Zeitskalen ist eine Näherung als linear aber noch vertretbar. Die Einschwingzeit scheint für die Berechnung nicht relevant zu sein, da in den ersten 100 Jahren kein deutlich größerer Unterschied des Ergebnisses zu den folgenden Zeitschritten zu erkennen ist.

Da unser berechneter Wert so stark vom Literaturwert abweicht, möchten wir die Parameter unserer Simulation anpassen um den Literaturwert zu erreichen. Dies verwenden wir für die folgenden Fits, um die Laufzeit zu verkürzen in dem wir die Berechnungen nur für die ersten 100 Jahre durchführen.

Shooting Methode und Variation der Flutbergmasse und Reibungskonstante¶

Die Reibung der Flutberge auf der Erdoberfläche sorgt für eine Verlangsamung der Erdrotation und einem größeren Orbit für den Mond. Wir haben für diese Berechnung zwei freie Parameter, die wir im folgenden anpassen wollen:

  • Die Masse der Flutberge $m_1$ und $m_2$ (die Masse des Ozeans)
  • Die Reibungskonstante $k$

Shooting Methode¶

Mit Hilfe der Shooting Methode kann man Randwertprobleme numerisch auf Anfangswertprobleme zurückführen. Es wird zu Anfang ein Intervall geraten in dem die mögliche Lösung liegt. Für die Mitte des Intervalls wird das Problem numerisch gelöst und mit dem Randwert verglichen. Es wird nun mit Hilfe des Bisektionverfahrens das Intervall solange halbiert bis die Lösung ausreichend genau gefunden wurde.

Variation der Masse der Flutberge¶

Wir modifizieren die Massen der Flutberge so, dass sie für die aktuellen Zunahme der Tageslänge pro 100 Jahre mit dem oben angegebenen Literaturwert $\tau$ übereinstimmt (hier bleibt $k$ unverändert).

In der Differentialgleichung für die Erdrotationsbeschleunigung sehen wir, dass diese linear mit der Masse der Flutberge zusammenhängt. Die berechnete Erdrotationsbeschleunigung ist zu schnell, deshalb müssen wir die Masse der Flutberge verkleinern. Wir wählen deshalb für die Shooting Methode als Suchintervall Null bis zur aktuellen Ozeanmasse. Es wäre möglich die Massen der Flutberge einzeln zu variieren, wir haben uns aber entschieden sie für die folgenden Berechnungen gleich zu setzen.

Variation der Reibungskonstante¶

Die verwendete Ozeanmasse entspricht bereits den reelen Begebenheiten, deshalb ist es nicht sinnvoll diese zu fitten. Alternativ können wir die Reibungskonstante verändern. Auch sie ist aus den gleichen Gründen zu groß und wir wählen als Suchintervall Null bis zum aktuellen Wert. Die Masse der Flutberge lassen wir hier unverändert.

Zum Ausführen der Shooting-Methode müssen wir noch eine Toleranz angeben. Der Literaturwert $\tau = 0.0021$ ist mit einer Genauigkeit von $10^{-5}$ angegeben, deshalb wählen wir diese als Toleranz.

In [ ]:
def rotation_perdiod_change_mOcean(mOcean: float):
    '''Calculates the change in rotation_period after 100 years, given a value for mOcean'''
    initial_values = iv_stable_orbit_4body_friction()
    # Change the mass of the two tides to each be half the value of mOcean for the simulation
    initial_values[-1][-1], initial_values[-1][-2] = mOcean/2, mOcean/2
    _, _, _, _, _, _, _, _, _, _, omega_E = four_body_problem_friction(*initial_values, k=2*10**(-12), t_max=100*365*24*3600) # Simulate 100 years
    return (2*pi/omega_E[-1]) - (2*pi/omega_E[0]) # Calculate the change in rotation period as described above

def rotation_period_change_k(k: float):
    '''Calculates the change in rotation_period after 100 years, given a value for k'''
    _, _, _, _, _, _, _, _, _, _, omega_E = four_body_problem_friction(*iv_stable_orbit_4body_friction(), k=k, t_max=100*365*24*3600) # Simulate 100 years
    return (2*pi/omega_E[-1]) - (2*pi/omega_E[0]) # Calculate the change in rotation period as described above

def shooting_method(func, target: float, val_min: float, val_max: float, tol: float):
    '''Find the value of x for which func(x) is in the interval [target-tol, target+tol]
    Search in the interval [val_min, val_max]'''

    # Check if the solution is in the passed interval
    if func(val_min) > target or func(val_max) < target:
        raise ValueError('The target is not in the passed interval')
    print("Target in interval")

    # Iterate until the solution is within the tolerance
    while True:
        x = (val_min + val_max)/2 # Calculate the midpoint of the interval
        func_val = func(x) # Calculate the value of the function at the midpoint
        if np.abs(func_val- target) < tol: # Solution found if the absolute distance to the solution is smaller than the tolerance
            return x
        elif func_val < target: # The lower bound can be raised
            val_min = x
        else: # The upper bound can be lowered
            val_max = x

# Fit the value of mOceans so that the change in rotation period is tau
if not fastExecution:
    mOzean_fit = shooting_method(func=rotation_perdiod_change_mOcean, target=tau, val_min=0, val_max=mOzean, tol=0.00001)
else:
    mOzean_fit = 2.13623046875e+16

# Fit the value of k so that the change in rotation period is tau
if not fastExecution:
    k_fit = shooting_method(func=rotation_period_change_k, target=tau, val_min=0, val_max=2*10**(-12), tol=0.00001)
else:
    k_fit = 3.06640625e-17

print(f'mOzean_fit: {mOzean_fit:.3g} kg')
print(f'k_fit: {k_fit:.3g} 1/m')
mOzean_fit: 2.14e+16 kg
k_fit: 3.07e-17 1/m

Ergebnisse Fit durch Shooting Methode¶

Wir erhalten wie erwartet durch den Fit für Ozeanmasse und Reibungskonstante kleinere Werte als unsere Startwerte.

Der gefittete Wert für die Ozeanmasse liegt aber mit $2.14 \cdot 10^{16} \mathrm{kg}$ deutlich unter dem vorher verwendeteten Literaturwert von $ 1.4 \cdot 10^{21} \mathrm{kg} $ Quelle: Ocean#Weigth.

Die Abweichungen lassen sich durch mehrere Näherungen in der Simulation erklären:

  • Korrektheit der Reibungskonstante
  • Die Sonne wird nicht berücksichtigt. Sie trägt einen kleinen Teil zu den Gezeiten bei.
  • Die Erde verformt sich durch die Gezeitenkräfte
  • Landmassen und ihre Effekte auf die Gezeiten werden nicht berücksichtigt

Nach den obigen Betrachtungen scheint es also sinnvoller statt der Ozeanmasse die Reibungskonstante zu fitten. Diese verhält sich auch viel mehr wie ein freier Parameter, da sie viele in der Realität auftretende Effekte zusammenfasst. (z.B. Landmassen, verschiedene Meerestiefen, Strömungen)

Vergleich der gefundenen Reibungskonstante mit der gegebenen¶

In [ ]:
if not fastExecution:
    # 4 body with friction  for 100 years: 11min
    t_kfit, x_E_kfit, y_E_kfit, x_M_kfit, y_M_kfit, phi_1_kfit, phi_2_kfit, phi_E_kfit, omega_1_kfit, omega_2_kfit, omega_E_kfit = four_body_problem_friction(*iv_stable_orbit_4body_friction(), k=k_fit, t_max=100*365*24*3600, rtol=1e-12, atol=1e-12)
    x_1_kfit, y_1_kfit, x_2_kfit, y_2_kfit = x_E_kfit + RErde*np.cos(phi_1_kfit), y_E_kfit + RErde*np.sin(phi_1_kfit), x_E_kfit + RErde*np.cos(phi_2_kfit), y_E_kfit + RErde*np.sin(phi_2_kfit)
In [ ]:
if not fastExecution:
    fig, axes = plt.subplots(nrows=2, ncols=3, figsize=(15, 10))
    plt.rc('axes', axisbelow=True)
    plt.suptitle("Verlgeich 4 Körper Problem $k=2*10^{-12}$ / $k=3*10^{-17}$")

    def plot_subplot(ax, lines, title, xlabel, ylabel):
        for x, y, label, color in lines:
            ax.plot(x, y, color, linewidth=1, label=label)
        ax.set_xlabel(xlabel)
        ax.set_ylabel(ylabel)
        ax.set_title(title)

        ax.grid('minor', 'minor', linestyle='-', linewidth=0.2)
        ax.grid('major', 'major', linestyle='-', linewidth=0.8)
        ax.minorticks_on()
        ax.legend()


    # Distance from the origin for earth and moon
    r_Earth_fl = np.sqrt(x_E_fl**2 + y_E_fl**2)
    r_Moon_fl = np.sqrt(x_M_fl**2 + y_M_fl**2)
    r_Earth_kfit = np.sqrt(x_E_kfit**2 + y_E_kfit**2)
    r_Moon_kfit = np.sqrt(x_M_kfit**2 + y_M_kfit**2)

    plot_subplot(axes[0][0], [(x_E_fl*10, y_E_fl*10, 'Erde (*10)', 'k'), (x_M_fl, y_M_fl, 'Mond', 'r'), (x_1_fl*10, y_1_fl*10, 'Flutberg 1 (*10)', 'b'), (x_2_fl*10, y_2_fl*10, 'Flutberg 2 (*10)', 'b')], "Bahnen des 4-Körper Problems $k=2*10^{-12}$", "x / m", "y / m")
    plot_subplot(axes[0][1], [(t_fl, r_Moon_fl, 'Mond', 'r')], "Abstand des Mondes zum Ursprung $k=2*10^{-12}$", "t / s", "r / m")
    plot_subplot(axes[0][2], [(t_fl, r_Earth_fl, 'Erde', 'k')], "Abstand der Erde zum Ursprung $k=2*10^{-12}$", "t / s", "r / m")
    plot_subplot(axes[1][0], [(x_E_kfit*10, y_E_kfit*10, 'Erde (*10)', 'k'), (x_M_kfit, y_M_kfit, 'Mond', 'r'), (x_1_kfit*10, y_1_kfit*10, 'Flutberg 1 (*10)', 'b'), (x_2_kfit*10, y_2_kfit*10, 'Flutberg 2 (*10)', 'b')], "Bahnen des 4-Körper Problems $k=3*10^{-17}$", "x / m", "y / m")
    plot_subplot(axes[1][1], [(t_kfit, r_Moon_kfit, 'Mond', 'r')], "Abstand des Mondes zum Ursprung $k=3*10^{-17}$", "t / s", "r / m")
    plot_subplot(axes[1][2], [(t_kfit, r_Earth_kfit, 'Erde', 'k')], "Abstand der Erde zum Ursprung $k=3*10^{-17}$", "t / s", "r / m")
    plt.tight_layout()
    plt.savefig("4body_friction_comp1_k_fit.png", dpi=300)
    plt.close()

    fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(10, 5))
    plt.rc('axes', axisbelow=True)
    plt.suptitle("Verlgeich 4 Körper Problem $k=2*10^{-12}$ / $k=3*10^{-17}$")

    plot_subplot(axes[0][0], [(t_fl, omega_1_fl, '$\omega_1$ Flutberg 1', 'k'), (t_fl, omega_2_fl, '$\omega_2$ Flutberg 2', 'r'), (t_fl, [2*pi/TMondBahn]*np.shape(t_fl)[0], 'Startwinkelgeschwindigkeit Flutberg', 'g')], "Winkelgeschwindigkeiten der Flutberge $k=2*10^{-12}$", "t / s", "$\omega$ / rad/s")
    plot_subplot(axes[0][1], [(t_fl, omega_E_fl, '$\omega_E$ Winkelgeschwindigkeit Erde', 'k'), (t_fl, [2*pi/TErdRotation]*np.shape(t_fl)[0], 'Startwinkelgeschwindigkeit Erde', 'g')], "Winkelgeschwindigkeit der Erde $k=2*10^{-12}$", "t / s", "$\omega$ / rad/s")
    plot_subplot(axes[1][0], [(t_kfit, omega_1_kfit, '$\omega_1$ Flutberg 1', 'k'), (t_kfit, omega_2_kfit, '$\omega_2$ Flutberg 2', 'r'), (t_kfit, [2*pi/TMondBahn]*np.shape(t_kfit)[0], 'Startwinkelgeschwindigkeit Flutberg', 'g')], "Winkelgeschwindigkeiten der Flutberge $k=3*10^{-17}$", "t / s", "$\omega$ / rad/s")
    plot_subplot(axes[1][1], [(t_kfit, omega_E_kfit, '$\omega_E$ Winkelgeschwindigkeit Erde', 'k'), (t_kfit, [2*pi/TErdRotation]*np.shape(t_kfit)[0], 'Startwinkelgeschwindigkeit Erde', 'g')], "Winkelgeschwindigkeit der Erde $k=3*10^{-17}$", "t / s", "$\omega$ / rad/s")

    plt.tight_layout()
    plt.savefig("4body_friction_comp2_k_fit.png", dpi=300)
    plt.close()

# Display image set the size to be 500x500 pixels
display(Image(filename='4body_friction_comp1_k_fit.png', width=1000, height=500))
display(Image(filename='4body_friction_comp2_k_fit.png', width=1000, height=500))

Da die gefittete Reibungskonstante um fünf Größenordnungen kleiner ist, sind die zuvor diskutierten Reibungseffekte sehr viel kleiner. Die Tageslängenänderung lässt sich jedoch trotzdem erkennen. Die anderen Effekte (Entfernung des Mondes von der Erde, Einschwingen), sind jedoch so klein, dass sie in den Plots nicht zu sehen sind.

Vielteilchenproblem mit Randbedingungen, intrinsischer Rotation und Reibung¶

Nun betrachten wir das Erde-Mond System mit $N$ Massepunkten auf der Erdoberfläche, welche die Wassermassen darstellen. Ziel ist es das Ausbilden der Flutberge zu simulieren. Für die Mond- und Erdkoordinaten verwenden wir die Differentialgleichung des Zweikörperproblems, als Näherung um den Rechenaufwand zu reduzieren. Die Näherung wird später genauer diskutiert. Die Kraft auf die Masspunkte ist gleich zur Kraft auf die zwei Flutberge mit Reibung, welche oben erläutert wurde. Die Differentialgleichung für $\omega_{\mathrm E}$ wird mit $N$ Massepunkten, analog zum Fall mit zwei Flutbergen, zu: $\newline$ \begin{equation*} \dot{\omega}_{\mathrm E} = \frac{5kR_{\mathrm E}}{2m_{\mathrm E}} \sum_{i = 1}^{N} m_i\lvert \omega_i - \omega_{\mathrm E} \rvert(\omega_i - \omega_{\mathrm E}) \end{equation*}

In [ ]:
N = 25 # Amount of tides

Wir verwenden wenn möglich die bereits implementierten Gleichungen. Die Winkel und Winkelbeschleunigungen der Flutberge werden als Liste übergeben, um $N$ leicht änden zu können.

In [ ]:
def eq_motion_Nbody(t, state, mass, k):
    '''This function calculates the derivatives of the state vector the 2 body problem with N Tides to be passed to solve_ivp. \\
    state: state vector has the following structure: [x_E, y_E, x_M, y_M, phi_E, vx_E, vy_E, vx_M, vy_M, omega_E, phi_1, phi_2, ..., phi_N, omega_1, omega_2, ..., omega_N]
    mass: list of masses e.g. [m_E, m_M, m1, m2, ..., mN]'''

    # Unpack the state vector and masses:
    x_E, y_E, x_M, y_M, phi_E, vx_E, vy_E, vx_M, vy_M, omega_E, phi_tide, omega_tide = *state[0:10], [state[10+i] for i in range(N)], [state[10+N+i] for i in range(N)]
    m_E, m_M, m_tide  = *mass[0:2], [mass[2+i] for i in range(N)]

    # Earth and moon acceleration from the 2 body problem:    
    _, _, _, _, ax_E, ay_E, ax_M, ay_M = eq_motion_2body(t, [x_E, y_E, x_M, y_M, vx_E, vy_E, vx_M, vy_M], [m_E, m_M])

    # Angular momentum of the earth:
    alpha_E = (5*k*RErde)/(2*m_E) * np.sum([m_tide[i]*abs(omega_tide[i] - omega_E)*(omega_tide[i] - omega_E) for i in range(N)])

    # The tides' cartesian coordinates:
    x_tide, y_tide = [x_E + RErde*np.cos(phi_tide[i]) for i in range(N)], [y_E + RErde*np.sin(phi_tide[i]) for i in range(N)]

    # The earth's angular acceleration:
    alpha_tide = [tides_acceleration_friction(x_M, y_M, x_tide[i], y_tide[i], phi_tide[i], ax_E, ay_E, omega_tide[i], omega_E, m_M, k) for i in range(N)]

    return [vx_E, vy_E, vx_M, vy_M, omega_E, ax_E, ay_E, ax_M, ay_M, alpha_E] + omega_tide + alpha_tide

Da für die Erd- und Mondkoordinaten die Zweikörper-Differentialgleichug verwendet wird wählen wir für diese dieselben Anfangsbedingungen. Die Massepunkte starten gleichmäßig verteilt auf der Erdoberfläche mit der mittleren Winkelgeschwindigkeit der Flutberge mit Reibung.

In [ ]:
if not fastExecution:
    _, _, _, _, _, _, _, _, omega_Niv, _, _ = four_body_problem_friction(*iv_stable_orbit_4body_friction(), k=k_fit, t_max=100*365*24*3600, rtol=1e-12, atol=1e-12)
    omega_i0 = np.mean(omega_Niv)
else:
    omega_i0 = 2.7107783823217206e-06
In [ ]:
def iv_stable_orbit_Nbody():
    '''Initial conditions for the Earth-Moon-NTides system in a stable orbit'''
    # Each of the tides gets an nth of the oceans mass 
    mass = [mErde, mMond] + [mOzean/N for i in range(N)]
    # Because the DEQ for the 2 body problem is used for the earth and moon movement, we can reuse the initial conditions (Tide mass evenly distributed around the earth)
    x0_Erde, x0_Mond, v0_Erde, v0_Mond, _ = iv_stable_orbit_2body()
    _, _, _, phi0_E, _, _, _, omega0_E, _ = iv_stable_orbit_4body_friction()

    phi0 = [(2*pi*i)/N for i in range(N)] # The tides all start spaced evenly around the earth
    #omega0 = [2*pi/TMondBahn]*N # The tides all start with the same angular velocity as the moon
    omega0 = [omega_i0]*N
    
    return [x0_Erde, x0_Mond, phi0_E, v0_Erde, v0_Mond, omega0_E, phi0, omega0, mass]

def N_body_problem(pos_body_1: list, pos_body_2: list, phi_E: float, vel_body_1: list, vel_body_2: list, omega_E: float, phi: list, omega: list, mass: list, k: float, t_max: float, t_start=0):
    solution = solve_ivp(fun=eq_motion_Nbody, t_span=[t_start, t_max], y0=[*pos_body_1, *pos_body_2, phi_E, *vel_body_1, *vel_body_2, omega_E, *phi, *omega], args=(mass, k), rtol=1e-12, atol=1e-12)
    x_E, y_E, x_M, y_M, phi_E, vx_E, vy_E, vx_M, vy_M, omega_E = solution.y[0:10]
    phi_tide = [solution.y[10+i] for i in range(N)]
    omega_tide = [solution.y[10+N+i] for i in range(N)]
    return [solution.t, x_E, y_E, x_M, y_M, phi_tide, omega_E, omega_tide]
In [ ]:
t, x_E, y_E, x_M, y_M, phi_tide, omega_E, omega_tide = N_body_problem(*iv_stable_orbit_Nbody(), k=k_fit, t_max=20*TMondBahn)
In [ ]:
# Cartesian coordinates of the tides:
x_tide, y_tide = [x_E + RErde*np.cos(phi_tide[i]) for i in range(N)], [y_E + RErde*np.sin(phi_tide[i]) for i in range(N)]

# Plot the movement of the earth, moon and tides
plt.rc('axes', axisbelow=True)
plt.figure(figsize=(5, 5))
plt.title("Bahnen des N-Flutberg-System")
plt.grid('minor', 'minor', linestyle='-', linewidth=0.2)
plt.grid('major', 'major', linestyle='-', linewidth=0.8)
plt.minorticks_on()
plt.xlabel("x")
plt.ylabel("y")
plt.plot(x_M/10, y_M/10, 'r', linewidth=1, label='Mond (/10)')
for i in range(N):
    plt.plot(x_tide[i], y_tide[i], 'b', linewidth=0.5, alpha=0.2)
plt.plot(x_E, y_E, 'k', linewidth=1, label='Erde')
plt.legend()
plt.show()
plt.close()

# Plot the angles of the tides over time, subtracting the angle of the moon
plt.rc('axes', axisbelow=True)
plt.figure(figsize=(5, 5))
plt.title("Winkel der Flutberge über Zeit, abzüglich $t\cdot \omega_0$")
plt.grid('minor', 'minor', linestyle='-', linewidth=0.2)
plt.grid('major', 'major', linestyle='-', linewidth=0.8)
plt.minorticks_on()
plt.xlabel("$t$ / s")
plt.ylabel("$\phi - t\cdot \omega_0$ / rad")

for i in range(N):
    plt.plot(t, phi_tide[i] - t*2*pi/TMondBahn, 'k', linewidth=0.7, alpha=0.7)

plt.show()
plt.close()

# Plot the angular velocity of the earth over time
plt.rc('axes', axisbelow=True)
plt.figure(figsize=(5, 5))
plt.title("Winkelgeschwindigkeit der Erde über Zeit")
plt.grid('minor', 'minor', linestyle='-', linewidth=0.2)
plt.grid('major', 'major', linestyle='-', linewidth=0.8)
plt.minorticks_on()
plt.xlabel("$t$ / s")
plt.ylabel("$\omega_E$ / rad")
plt.plot(t, omega_E, 'k', linewidth=1, label='$\omega_E$')

plt.show()
plt.close()

# Animation of the N body problem
if not fastExecution:
    fig = plt.figure(figsize=(5, 5))
    ax = plt.axes(xlim=(-4e8, 4e8), ylim=(-4e8, 4e8))
    ax.set_title("Animation des N-Flutberg-Systems")
    earth_line, = ax.plot([], [], marker='o', lw=0.1, color='blue', label='Erde (*10)')
    moon_line, = ax.plot([], [], marker='o', lw=0.1, color='red', label='Mond')
    tide_lines = [ax.plot([], [], marker='o', lw=0.1, color='black', label='Flutberg {i} (*10)', alpha=0.15)[0] for i in range(N)]

    def animate(i):
        earth_line.set_data([x_E[i]*10], [y_E[i]*10])
        moon_line.set_data([x_M[i]], [y_M[i]])
        for j in range(N):
            tide_lines[j].set_data([x_tide[j][i]*10], [y_tide[j][i]*10])
        return earth_line, moon_line, *tide_lines

    anim = FuncAnimation(fig, animate, init_func=None, frames=t.shape[0], interval=30, blit=True)
    anim.save('Erde_Mond_NFlutberge.gif', writer='pillow')
    plt.close()

display(Image(data=open('Erde_Mond_NFlutberge.gif','rb').read(), format='png'))

Die Flutberge bilden sich aus, sind aber nicht sehr stabil. Ein paar Massen sind immer außerhalb der Flutberge, die Mehrzahl sammelt sich jedoch dort, wo die Flutberge sein sollten. Sie oszillieren auch unterschiedlich stark um die theoretische Position der Flutberge (abhängig vom Startpunkt).

Die Näherung ist gerechtfertigt, da einerseits die Gesammtmasse der Ozeane relativ zur Erd- und Mondmasse recht klein ist. Andererseits verkleinert sich der Effekt der Flutberge auf die Erd- und Mondkoordinaten durch die Verteilung der Masse. Auf die Erde wirken die Flutberge so, dass ihre Radialbeschleunigung gegen die Erde "drückt" beziehungsweise an ihr "zieht". Da die Massen hier verteilt sind mitteln sich diese Beschleunigungen zum Teil, daher kann dieser Effekt besser vernachlässigt werden als im Fall von 2 Flutbergen. Auf den Mond wirkt die Anziehungskraft der Flutmassen, welche bewirkt, dass der Mond ein Drehmoment erfährt. Dieser Effekt ist relativ klein, wodurch der Einfluss auf die Bewegung der Flutberge ebenfalls sehr klein ist. Dadurch entfernt sich der Mond hier nicht von der Erde, was für große Zeiten ein relevanter Effekt ist.

In [ ]:
# Two tides to compare to the previous result:
if not fastExecution:
    N = 2
    t_N, x_E_N, y_E_N, x_M_N, y_M_N, phi_tides_N, omega_E_N, omega_tides_N = N_body_problem(*iv_stable_orbit_Nbody(), k=k_fit, t_max=100*365*24*3600)
    phi_1_N, phi_2_N = phi_tides_N[0], phi_tides_N[1]
    omega_1_N, omega_2_N = omega_tides_N[0], omega_tides_N[1]
    x_1_N, y_1_N, x_2_N, y_2_N = x_E_N + RErde*np.cos(phi_1_N), y_E_N + RErde*np.sin(phi_1_N), x_E_N + RErde*np.cos(phi_2_N), y_E_N + RErde*np.sin(phi_2_N)
In [ ]:
if not fastExecution:
    fig, axes = plt.subplots(nrows=2, ncols=3, figsize=(15, 10))
    plt.rc('axes', axisbelow=True)
    plt.suptitle("Verlgeich 4 Körper Problem Exakte DGL / N=2 Näherung")

    def plot_subplot(ax, lines, title, xlabel, ylabel):
        for x, y, label, color in lines:
            ax.plot(x, y, color, linewidth=1, label=label)
        ax.set_xlabel(xlabel)
        ax.set_ylabel(ylabel)
        ax.set_title(title)

        ax.grid('minor', 'minor', linestyle='-', linewidth=0.2)
        ax.grid('major', 'major', linestyle='-', linewidth=0.8)
        ax.minorticks_on()
        ax.legend()


    # Distance from the origin for earth and moon
    r_Earth_kfit = np.sqrt(x_E_kfit**2 + y_E_kfit**2)
    r_Moon_kfit = np.sqrt(x_M_kfit**2 + y_M_kfit**2)
    r_Earth_N = np.sqrt(x_E_N**2 + y_E_N**2)
    r_Moon_N = np.sqrt(x_M_N**2 + y_M_N**2)

    plot_subplot(axes[0][0], [(x_E_kfit*10, y_E_kfit*10, 'Erde (*10)', 'k'), (x_M_kfit, y_M_kfit, 'Mond', 'r'), (x_1_kfit*10, y_1_kfit*10, 'Flutberg 1 (*10)', 'b'), (x_2_kfit*10, y_2_kfit*10, 'Flutberg 2 (*10)', 'b')], "Bahnen des 4-Körper Problems Exakte DGL", "x / m", "y / m")
    plot_subplot(axes[0][1], [(t_kfit, r_Moon_kfit, 'Mond', 'r')], "Abstand des Mondes zum Ursprung Exakte DGL", "t / s", "r / m")
    plot_subplot(axes[0][2], [(t_kfit, r_Earth_kfit, 'Erde', 'k')], "Abstand der Erde zum Ursprung Exakte DGL", "t / s", "r / m")
    plot_subplot(axes[1][0], [(x_E_N*10, y_E_N*10, 'Erde (*10)', 'k'), (x_M_N, y_M_N, 'Mond', 'r'), (x_1_N*10, y_1_N*10, 'Flutberg 1 (*10)', 'b'), (x_2_N*10, y_2_N*10, 'Flutberg 2 (*10)', 'b')], "Bahnen des 4-Körper Problems N=2 Näherung", "x / m", "y / m")
    plot_subplot(axes[1][1], [(t_N, r_Moon_N, 'Mond', 'r')], "Abstand des Mondes zum Ursprung N=2 Näherung", "t / s", "r / m")
    plot_subplot(axes[1][2], [(t_N, r_Earth_N, 'Erde', 'k')], "Abstand der Erde zum Ursprung N=2 Näherung", "t / s", "r / m")
    plt.tight_layout()
    plt.savefig("4body_friction_comp1_N=2.png", dpi=300)
    plt.close()

    fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(10, 5))
    plt.rc('axes', axisbelow=True)
    plt.suptitle("Verlgeich 4 Körper Problem Exakte DGL / N=2 Näherung")

    plot_subplot(axes[0][0], [(t_kfit, omega_1_kfit, '$\omega_1$ Flutberg 1', 'k'), (t_kfit, omega_2_kfit, '$\omega_2$ Flutberg 2', 'r'), (t_kfit, [2*pi/TMondBahn]*np.shape(t_kfit)[0], 'Startwinkelgeschwindigkeit Flutberg', 'g')], "Winkelgeschwindigkeiten der Flutberge Exakte DGL", "t / s", "$\omega$ / rad/s")
    plot_subplot(axes[0][1], [(t_kfit, omega_E_kfit, '$\omega_E$ Winkelgeschwindigkeit Erde', 'k'), (t_kfit, [2*pi/TErdRotation]*np.shape(t_kfit)[0], 'Startwinkelgeschwindigkeit Erde', 'g')], "Winkelgeschwindigkeit der Erde Exakte DGL", "t / s", "$\omega$ / rad/s")
    plot_subplot(axes[1][0], [(t_N, omega_1_N, '$\omega_1$ Flutberg 1', 'k'), (t_N, omega_2_N, '$\omega_2$ Flutberg 2', 'r'), (t_N, [omega_i0]*np.shape(t_N)[0], 'Startwinkelgeschwindigkeit Flutberg', 'g')], "Winkelgeschwindigkeiten der Flutberge N=2 Näherung", "t / s", "$\omega$ / rad/s")
    plot_subplot(axes[1][1], [(t_N, omega_E_N, '$\omega_E$ Winkelgeschwindigkeit Erde', 'k'), (t_N, [2*pi/TErdRotation]*np.shape(t_N)[0], 'Startwinkelgeschwindigkeit Erde', 'g')], "Winkelgeschwindigkeit der Erde N=2 Näherung", "t / s", "$\omega$ / rad/s")

    plt.tight_layout()
    plt.savefig("4body_friction_comp2_N=2.png", dpi=300)
    plt.close()

# Display image set the size to be 500x500 pixels
display(Image(filename='4body_friction_comp1_N=2.png', width=1000, height=500))
display(Image(filename='4body_friction_comp2_N=2.png', width=1000, height=500))
In [ ]:
#Show the incerase of the moons orbit
if not fastExecution:
    #Moon as seen from Earth:
    x = x_M_N - x_E_N
    y = y_M_N - y_E_N

    fig = plt.figure(figsize=(12,6))
    fig.suptitle('Veränderung der Mondbahn über 100 Jahre')
    ax = fig.add_subplot(121)
    ax2 = fig.add_subplot(122)

    ax.set_title('Mond im Erdsystem')
    ax.set_xlabel("$x$ / m")
    ax.set_ylabel("$y$ / m")
    ax.set_aspect('equal')
    ax.plot(x, y, color='r', linewidth=0.5, label='Mond')

    rect = plt.Rectangle((3.75e8, -1e5), 1.5e7, 1e7, facecolor="black", alpha=0.7, label='Ausschnitt')
    ax.add_patch(rect)
    ax.legend()

    ax2.set_title('Ausschnitt')
    ax2.set_xlabel("$x$ / m")
    ax2.set_ylabel("$y$ / m")
    ax2.set_aspect('equal')
    ax2.plot(x, y, 'r', linewidth=0.1, label='Mond')
    ax2.scatter(x[0], y[0], c='k', s=10, label='Startpunkt')
    ax2.set_xlim([3.8358e8, 3.8391e8])
    ax2.set_ylim([-1e5, 1e5])
    ax2.legend()

    plt.tight_layout()
    plt.savefig("4body_friction_Moon_N=2.png", dpi=300)
    plt.close()

display(Image(filename='4body_friction_Moon_N=2.png', width=1000, height=2000))
In [ ]:
if not fastExecution:
    N = 2
    t, x_E, y_E, x_M, y_M, phi_tide, omega_E, omega_tide = N_body_problem(*iv_stable_orbit_Nbody(), k=k_fit, t_max=5*TMondBahn)
In [ ]:
# Cartesian coordinates of the tides:
x_tide, y_tide = [x_E + RErde*np.cos(phi_tide[i]) for i in range(N)], [y_E + RErde*np.sin(phi_tide[i]) for i in range(N)]

# Animation of the N body problem
if not fastExecution: # Don't render if not necessary
    fig = plt.figure(figsize=(5, 5)) # Square figure, so the circles are actually circles
    ax = plt.axes(xlim=(-4e8, 4e8), ylim=(-4e8, 4e8))
    earth_line, = ax.plot([], [], marker='o', lw=0.1, color='blue', label='Erde (*10)')
    moon_line, = ax.plot([], [], marker='o', lw=0.1, color='red', label='Mond')
    tide_lines = [ax.plot([], [], marker='o', lw=0.1, color='black', label='Flutberg {i} (*10)', alpha=0.15)[0] for i in range(N)]
    # ax.legend()

    def animate(i):
        earth_line.set_data([x_E[i]*10], [y_E[i]*10])
        moon_line.set_data([x_M[i]], [y_M[i]])
        for j in range(N):
            tide_lines[j].set_data([x_tide[j][i]*10], [y_tide[j][i]*10])
        return earth_line, moon_line, *tide_lines

    anim = FuncAnimation(fig, animate, init_func=None, frames=t.shape[0], interval=30, blit=True)
    anim.save('Erde_Mond_N=2Flutberge.gif', writer='pillow')
    plt.close()

display(Image(data=open('Erde_Mond_N=2Flutberge.gif','rb').read(), format='png'))
In [ ]:
if not fastExecution:
    #30min
    _, _,_, _, _, _, omega_E, _ = N_body_problem(*iv_stable_orbit_Nbody(), k=k_fit, t_max=100*365*24*3600)

    # The change in the angular velocity
    angular_velocity_change = (omega_E[0] - omega_E[-1])

    # w = 2*pi/T -> T = 2*pi/w
    rotation_period_at_start = (2*pi/omega_E[0])
    rotation_period_at_end = (2*pi/omega_E[-1])
    print(f'rotation_period_at_start: {rotation_period_at_start}, rotation_period_at_end: {rotation_period_at_end}')
    print(f'rotation_period_at_start [h]: {rotation_period_at_start/3600}, rotation_period_at_end [h]: {rotation_period_at_end/3600}')

    print(f'tau_lit: {tau}, tau_calc: {rotation_period_at_end - rotation_period_at_start}')
else:
    print(f'tau_lit: {tau}, tau_calc: {0.0021018376137362793}')
tau_lit: 0.0021, tau_calc: 0.0021018376137362793

Die größten Unterschiede zur exakten Differentialgleichung sind:
Der Mond entfernt sich nicht von der Erde, da die Gravitation der Flutberge auf den Mond vernachlässigt wird. Auch die Schwerpunktsbewegung unterscheidet sich, was aufgrund der entsprechend gewählten Anfangsbedingungen nicht zu sehen ist.

Die Bewegung der Flutberge bleibt sehr ähnlich, um diese zu betrachten ist die Näherung also gut. Daher bleibt auch die Verlangsamung der Erdrotation fast gleich. Die Winkelgeschwindigkeiten unterscheiden sich nur so stark, weil wir unterschiedliche Anfangswerte für sie gewählt haben. Insgesamt sehen sich die Lösungen sehr ähnlich, und solange man nicht an der Änderung der Mondbahn interessiet ist, ist die Näherung auch gerechtfertigt. Für sehr lange Zeiten ist sie dementsprechend nicht ausreichend.